R Notebook

library(tidyverse)

── Attaching packages ───────────────────────────────── tidyverse 1.2.1 ──
✔ ggplot2 3.2.1     ✔ readr   1.3.1
✔ tibble  2.1.3     ✔ purrr   0.3.2
✔ tidyr   1.0.0     ✔ stringr 1.4.0
✔ ggplot2 3.2.1     ✔ forcats 0.4.0
── Conflicts ──────────────────────────────────── tidyverse_conflicts() ──
✖ tidyr::extract() masks raster::extract()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()
✖ dplyr::select()  masks raster::select()

library(sf)
library(raster)
library(spData)
library(spDataLarge)
library(lwgeom)

Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3

knitr::opts_chunk$set()

Vector Data

An introduction to simple features

names(world)

 [1] "iso_a2"    "name_long" "continent" "region_un" "subregion"
 [6] "type"      "area_km2"  "pop"       "lifeExp"   "gdpPercap"
[11] "geom"     

world

Simple feature collection with 177 features and 10 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

plot(world, max.plot = 10)

summary(world["lifeExp"])

    lifeExp                 geom    
 Min.   :50.62   MULTIPOLYGON :177  
 1st Qu.:64.96   epsg:4326    :  0  
 Median :72.87   +proj=long...:  0  
 Mean   :70.85                      
 3rd Qu.:76.78                      
 Max.   :83.59                      
 NA's   :10                         

world$geom[[1]]

MULTIPOLYGON (((180 -16.06713, 180 -16.55522, 179.3641 -16.80135, 178.7251 -17.01204, 178.5968 -16.63915, 179.0966 -16.43398, 179.4135 -16.37905, 180 -16.06713)), ((178.1256 -17.50481, 178.3736 -17.33992, 178.7181 -17.62846, 178.5527 -18.15059, 177.9327 -18.28799, 177.3815 -18.16432, 177.285 -17.72465, 177.6709 -17.38114, 178.1256 -17.50481)), ((-179.7933 -16.02088, -179.9174 -16.50178, -180 -16.55522, -180 -16.06713, -179.7933 -16.02088)))

world_mini = world[1:2, 1:3]
world_mini

Simple feature collection with 2 features and 3 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -18.28799 xmax: 180 ymax: -0.95
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

library(sp)
world_sp = as(world, Class = "Spatial")
world_sp

class       : SpatialPolygonsDataFrame 
features    : 177 
extent      : -180, 180, -90, 83.64513  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
variables   : 10
names       : iso_a2,   name_long,     continent,               region_un,      subregion,              type,         area_km2,        pop,          lifeExp,        gdpPercap 
min values  :     AE, Afghanistan,        Africa,                  Africa,     Antarctica,           Country, 2416.87048266498,      56295,           50.621, 597.135168986395 
max values  :     ZW,    Zimbabwe, South America, Seven seas (open ocean), Western Europe, Sovereign country, 17018507.4094666, 1364270000, 83.5878048780488,  120860.06755829 

world_sf = st_as_sf(world_sp)
world_sf

Simple feature collection with 177 features and 10 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
First 10 features:
   iso_a2        name_long     continent region_un          subregion
1      FJ             Fiji       Oceania   Oceania          Melanesia
2      TZ         Tanzania        Africa    Africa     Eastern Africa
3      EH   Western Sahara        Africa    Africa    Northern Africa
4      CA           Canada North America  Americas   Northern America
5      US    United States North America  Americas   Northern America
6      KZ       Kazakhstan          Asia      Asia       Central Asia
7      UZ       Uzbekistan          Asia      Asia       Central Asia
8      PG Papua New Guinea       Oceania   Oceania          Melanesia
9      ID        Indonesia          Asia      Asia South-Eastern Asia
10     AR        Argentina South America  Americas      South America
                type    area_km2       pop  lifeExp gdpPercap
1  Sovereign country    19289.97    885806 69.96000  8222.254
2  Sovereign country   932745.79  52234869 64.16300  2402.099
3      Indeterminate    96270.60        NA       NA        NA
4  Sovereign country 10036042.98  35535348 81.95305 43079.143
5            Country  9510743.74 318622525 78.84146 51921.985
6  Sovereign country  2729810.51  17288285 71.62000 23587.338
7  Sovereign country   461410.26  30757700 71.03900  5370.866
8  Sovereign country   464520.07   7755785 65.23000  3709.082
9  Sovereign country  1819251.33 255131116 68.85600 10003.089
10 Sovereign country  2784468.59  42981515 76.25200 18797.548
                         geometry
1  MULTIPOLYGON (((180 -16.067...
2  MULTIPOLYGON (((33.90371 -0...
3  MULTIPOLYGON (((-8.66559 27...
4  MULTIPOLYGON (((-122.84 49,...
5  MULTIPOLYGON (((-122.84 49,...
6  MULTIPOLYGON (((87.35997 49...
7  MULTIPOLYGON (((55.96819 41...
8  MULTIPOLYGON (((141.0002 -2...
9  MULTIPOLYGON (((141.0002 -2...
10 MULTIPOLYGON (((-68.63401 -...

Basic map making

plot(world[3:6])

plot(world["pop"])

world_asia <- world %>% 
  filter(continent == "Asia")
world_asia

Simple feature collection with 47 features and 10 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 26.04335 ymin: -10.35999 xmax: 145.5431 ymax: 55.38525
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

asia = st_union(world_asia)
asia

Geometry set for 1 feature 
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 26.04335 ymin: -10.35999 xmax: 145.5431 ymax: 55.38525
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

MULTIPOLYGON (((120.295 -10.25865, 118.9678 -9....

plot(world["pop"], reset = FALSE)
plot(asia, add = TRUE, col = "red")

plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)

st_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

plot(st_geometry(world_cents), add = TRUE, cex = cex)

india = world[world$name_long == "India", ]
plot(st_geometry(india), expandBB = c(0, 0.2, 0.1, 1), col = "gray", lwd = 3)
plot(world_asia[0], add = TRUE)

Geometry tyeps

Simple feature geometries (sfg)

st_point(c(5, 2))                 # XY point

POINT (5 2)

st_point(c(5, 2, 3))              # XYZ point

POINT Z (5 2 3)

st_point(c(5, 2, 1), dim = "XYM") # XYM point

POINT M (5 2 1)

st_point(c(5, 2, 3, 1))           # XYZM point

POINT ZM (5 2 3 1)

# the rbind function simplifies the creation of matrices

## MULTIPOINT
multipoint_matrix = rbind(c(5, 2), c(1, 3), c(3, 4), c(3, 2))
st_multipoint(multipoint_matrix)

MULTIPOINT (5 2, 1 3, 3 4, 3 2)

## LINESTRING
linestring_matrix = rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2))
st_linestring(linestring_matrix)

LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2)

## POLYGON
polygon_list = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
st_polygon(polygon_list)

POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))

## POLYGON with a hole
polygon_border = rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))
polygon_hole = rbind(c(2, 4), c(3, 4), c(3, 3), c(2, 3), c(2, 4))
polygon_with_hole_list = list(polygon_border, polygon_hole)
st_polygon(polygon_with_hole_list)

POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5), (2 4, 3 4, 3 3, 2 3, 2 4))

## MULTILINESTRING
multilinestring_list = list(
  rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
  rbind(c(1, 2), c(2, 4))
  )
st_multilinestring((multilinestring_list))

MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 2, 2 4))

## MULTIPOLYGON
multipolygon_list = list(
  list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5))),
  list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2)))
  )
st_multipolygon(multipolygon_list)

MULTIPOLYGON (((1 5, 2 2, 4 1, 4 4, 1 5)), ((0 2, 1 2, 1 3, 0 3, 0 2)))

## GEOMETRYCOLLECTION
gemetrycollection_list = list(st_multipoint(multipoint_matrix),
                              st_linestring(linestring_matrix))
st_geometrycollection(gemetrycollection_list)

GEOMETRYCOLLECTION (MULTIPOINT (5 2, 1 3, 3 4, 3 2), LINESTRING (1 5, 4 4, 4 1, 2 2, 3 2))

Simple feature columns (sfc)

# sfc POINT
point1 = st_point(c(5, 2))
point2 = st_point(c(1, 3))
points_sfc = st_sfc(point1, point2)
points_sfc

Geometry set for 2 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 2 xmax: 5 ymax: 3
epsg (SRID):    NA
proj4string:    NA

POINT (5 2)
POINT (1 3)

# sfc POLYGON
polygon_list1 = list(rbind(c(1, 5), c(2, 2), c(4, 1), c(4, 4), c(1, 5)))
polygon1 = st_polygon(polygon_list1)

polygon_list2 = list(rbind(c(0, 2), c(1, 2), c(1, 3), c(0, 3), c(0, 2)))
polygon2 = st_polygon(polygon_list2)

polygon_sfc = st_sfc(polygon1, polygon2)
polygon_sfc

Geometry set for 2 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 0 ymin: 1 xmax: 4 ymax: 5
epsg (SRID):    NA
proj4string:    NA

POLYGON ((1 5, 2 2, 4 1, 4 4, 1 5))
POLYGON ((0 2, 1 2, 1 3, 0 3, 0 2))

# sfc MULTILINESTRING
multilinestring_list1 = list(
  rbind(c(1, 5), c(4, 4), c(4, 1), c(2, 2), c(3, 2)), 
  rbind(c(1, 2), c(2, 4))
  )
multilinestring1 = st_multilinestring((multilinestring_list1))

multilinestring_list2 = list(
  rbind(c(2, 9), c(7, 9), c(5, 6), c(4, 7), c(2, 7)), 
  rbind(c(1, 7), c(3, 8))
  )
multilinestring2 = st_multilinestring((multilinestring_list2))

multilinestring_sfc = st_sfc(multilinestring1, multilinestring2)
multilinestring_sfc

Geometry set for 2 features 
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 1 ymin: 1 xmax: 7 ymax: 9
epsg (SRID):    NA
proj4string:    NA

MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 ...
MULTILINESTRING ((2 9, 7 9, 5 6, 4 7, 2 7), (1 ...

# sfc GEOMETRY
point_multilinestring_sfc = st_sfc(point1, multilinestring1)
point_multilinestring_sfc

Geometry set for 2 features 
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: 1 ymin: 1 xmax: 5 ymax: 5
epsg (SRID):    NA
proj4string:    NA

POINT (5 2)
MULTILINESTRING ((1 5, 4 4, 4 1, 2 2, 3 2), (1 ...

points_sfc

Geometry set for 2 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 2 xmax: 5 ymax: 3
epsg (SRID):    NA
proj4string:    NA

POINT (5 2)
POINT (1 3)

st_crs(points_sfc)

Coordinate Reference System: NA

# EPSG definition
st_sfc(point1, point2, crs = 4326)

Geometry set for 2 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 2 xmax: 5 ymax: 3
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

POINT (5 2)
POINT (1 3)

# PROJ4STRING definition
st_sfc(point1, point2, crs = "+proj=longlat +datum=WGS84 +no_defs")

Geometry set for 2 features 
geometry type:  POINT
dimension:      XY
bbox:           xmin: 1 ymin: 2 xmax: 5 ymax: 3
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

POINT (5 2)
POINT (1 3)

The sf class

lnd_point = st_point(c(0.1, 51.5))                 # sfg object
lnd_geom = st_sfc(lnd_point, crs = 4326)           # sfc object
lnd_attrib = data.frame(                           # data.frame object
  name = "London",
  temperature = 25,
  date = as.Date("2017-06-21")
  )
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object
lnd_sf

Simple feature collection with 1 feature and 3 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 0.1 ymin: 51.5 xmax: 0.1 ymax: 51.5
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
    name temperature       date         geometry
1 London          25 2017-06-21 POINT (0.1 51.5)

class(lnd_sf)

[1] "sf"         "data.frame"

Raster data

An introduction to raster

new_raster

class      : RasterLayer 
dimensions : 457, 465, 212505  (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333  (x, y)
extent     : -113.2396, -112.8521, 37.13208, 37.51292  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /Library/Frameworks/R.framework/Versions/3.6/Resources/library/spDataLarge/raster/srtm.tif 
names      : srtm 
values     : 1024, 2892  (min, max)

dim(new_raster)

[1] 457 465   1

ncell(new_raster)

[1] 212505

res(new_raster)

[1] 0.0008333333 0.0008333333

extent(new_raster)

class      : Extent 
xmin       : -113.2396 
xmax       : -112.8521 
ymin       : 37.13208 
ymax       : 37.51292 

crs(new_raster)

CRS arguments:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 

inMemory(new_raster)

[1] FALSE

Basic map making

Raster clsses

raster::writeFormats()

      name        long_name                                
 [1,] "raster"    "R-raster"                               
 [2,] "SAGA"      "SAGA GIS"                               
 [3,] "IDRISI"    "IDRISI"                                 
 [4,] "IDRISIold" "IDRISI (img/doc)"                       
 [5,] "BIL"       "Band by Line"                           
 [6,] "BSQ"       "Band Sequential"                        
 [7,] "BIP"       "Band by Pixel"                          
 [8,] "ascii"     "Arc ASCII"                              
 [9,] "CDF"       "NetCDF"                                 
[10,] "big"       "big.matrix"                             
[11,] "ADRG"      "ARC Digitized Raster Graphics"          
[12,] "BMP"       "MS Windows Device Independent Bitmap"   
[13,] "BT"        "VTP .bt (Binary Terrain) 1.3 Format"    
[14,] "BYN"       "Natural Resources Canada's Geoid"       
[15,] "CTable2"   "CTable2 Datum Grid Shift"               
[16,] "EHdr"      "ESRI .hdr Labelled"                     
[17,] "ELAS"      "ELAS"                                   
[18,] "ENVI"      "ENVI .hdr Labelled"                     
[19,] "ERS"       "ERMapper .ers Labelled"                 
[20,] "GPKG"      "GeoPackage"                             
[21,] "GS7BG"     "Golden Software 7 Binary Grid (.grd)"   
[22,] "GSBG"      "Golden Software Binary Grid (.grd)"     
[23,] "GTiff"     "GeoTIFF"                                
[24,] "GTX"       "NOAA Vertical Datum .GTX"               
[25,] "HFA"       "Erdas Imagine Images (.img)"            
[26,] "IDA"       "Image Data and Analysis"                
[27,] "ILWIS"     "ILWIS Raster Map"                       
[28,] "INGR"      "Intergraph Raster"                      
[29,] "ISCE"      "ISCE raster"                            
[30,] "ISIS2"     "USGS Astrogeology ISIS cube (Version 2)"
[31,] "ISIS3"     "USGS Astrogeology ISIS cube (Version 3)"
[32,] "KRO"       "KOLOR Raw"                              
[33,] "LAN"       "Erdas .LAN/.GIS"                        
[34,] "Leveller"  "Leveller heightfield"                   
[35,] "MBTiles"   "MBTiles"                                
[36,] "MRF"       "Meta Raster Format"                     
[37,] "netCDF"    "Network Common Data Format"             
[38,] "NGW"       "NextGIS Web"                            
[39,] "NITF"      "National Imagery Transmission Format"   
[40,] "NTv2"      "NTv2 Datum Grid Shift"                  
[41,] "NWT_GRD"   "Northwood Numeric Grid Format .grd/.tab"
[42,] "PAux"      "PCI .aux Labelled"                      
[43,] "PCIDSK"    "PCIDSK Database File"                   
[44,] "PCRaster"  "PCRaster Raster File"                   
[45,] "PDF"       "Geospatial PDF"                         
[46,] "PDS4"      "NASA Planetary Data System 4"           
[47,] "PNM"       "Portable Pixmap Format (netpbm)"        
[48,] "RMF"       "Raster Matrix Format"                   
[49,] "ROI_PAC"   "ROI_PAC raster"                         
[50,] "RRASTER"   "R Raster"                               
[51,] "RST"       "Idrisi Raster A.1"                      
[52,] "SAGA"      "SAGA GIS Binary Grid (.sdat, .sg-grd-z)"
[53,] "SGI"       "SGI Image File Format 1.0"              
[54,] "Terragen"  "Terragen heightfield"                   

r_brick

class      : RasterBrick 
dimensions : 1428, 1128, 1610784, 4  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : /Library/Frameworks/R.framework/Versions/3.6/Resources/library/spDataLarge/raster/landsat.tif 
names      : landsat.1, landsat.2, landsat.3, landsat.4 
min values :      7550,      6404,      5678,      5252 
max values :     19071,     22051,     25780,     31961 

nlayers(r_brick)

[1] 4

plot(r_brick)

raster_on_disk = raster(r_brick, layer = 1)
plot(raster_on_disk)

raster_in_memory = raster(xmn = 301905, xmx = 335745,
                          ymn = 4111245, ymx = 4154085, 
                          res = 30)
values(raster_in_memory) = sample(seq_len(ncell(raster_in_memory)))
crs(raster_in_memory) = crs(raster_on_disk)
plot(raster_in_memory)

r_stack = stack(raster_in_memory, raster_on_disk)
r_stack

class      : RasterStack 
dimensions : 1428, 1128, 1610784, 2  (nrow, ncol, ncell, nlayers)
resolution : 30, 30  (x, y)
extent     : 301905, 335745, 4111245, 4154085  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      :   layer, landsat.1 
min values :       1,      7550 
max values : 1610784,     19071 

Coordinate Reference System

Geographic coordinate system

Projected coordinate reference systems

CRS in R

vector_filepath = system.file("vector/zion.gpkg", package = "spDataLarge")
new_vector = st_read(vector_filepath)

Reading layer `zion' from data source `/Library/Frameworks/R.framework/Versions/3.6/Resources/library/spDataLarge/vector/zion.gpkg' using driver `GPKG'
Simple feature collection with 1 feature and 11 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 302903.1 ymin: 4112244 xmax: 334735.5 ymax: 4153087
epsg (SRID):    NA
proj4string:    +proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

st_crs(new_vector)

Coordinate Reference System:
  No EPSG code
  proj4string: "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

new_vector <- st_set_crs(new_vector, 4326)

st_crs<- : replacing crs does not reproject data; use st_transform for that

projection(new_raster) # get CRS

[1] "+proj=utm +zone=12 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

projection(new_raster) = "+ proj=utm + zone=12 + ellps=GRS80 + towgs84=0,0,0,0,0,0,0 + units=m + no_defs" # set CRS

+ proj=utm + zone=12 + ellps=GRS80 + towgs84=0,0,0,0,0,0,0 + units=m + no_defs is not a valid PROJ.4 CRS string

Units

st_area(luxembourg)

2416870483 [m^2]

st_area(luxembourg) / 1000000

2416.87 [m^2]

units::set_units(st_area(luxembourg), km^2)

2416.87 [km^2]

res(new_raster)

[1] 0.0008333333 0.0008333333

res(repr)

[1] 0.000833 0.000833